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ABSTRACT 

The recent development of a new minimum mass solar nebula, under the assumption that the 
giant planets formed in the compact configuration of the Nice model, has shed new light on planet 
formation in the solar system. Desch previously found that a steady state protoplanetary disk with an 
outer boundary truncated by photoevaporation by an external massive star would have a steep surface 
density profile. In a completely novel way, we have adapted numerical methods for solving propagating 
phase change problems to astrophysical disks. We find that a one-dimensional time-dependent disk 
model that self-consistently tracks the location of the outer boundary produces shallower profiles 
than those predicted for a steady state disk. The resulting surface density profiles have a radial 

dependence of S(r) oc r -° m with a power law exponent that in some models becomes as large 
as ~ S(r) oc r~ 2A . The evolutionary timescales of the model disks can be sped up or slowed down 
by altering the amount of far-ultraviolet flux or the viscosity parameter a. Slowing the evolutionary 
timescale by decreasing the incident far ultraviolet flux, or similarly by decreasing a, can help to grow 
planets more rapidly, but at the cost of decreased migration timescales. Although they similarly affect 
relevant timescales, changes in the far ultraviolet flux or a produce disks with drastically different 
outer radii. Despite their differences, these disks are all characterized by outward mass transport, 
mass loss at the outer edge, and a truncated outer boundary. The transport of mass from small to 
large radii can potentially prevent the rapid inward migration of Jupiter and Saturn, while at the same 
time supply enough mass to the outer regions of the disk for the formation of Uranus and Neptune. 
Subject headings: accretion, accretion disks - planet-disk interactions - planets and satellites: dynam- 
ical evolution and stability - planets and satellites: formation - protoplanetary 
disks 



1. INTRODUCTION 

The development of the minimum mass solar nebula 
(MMSN) was a critical step in understanding the ori- 
gin of the solar system . The MMSN w as developed b y 
iWeidenschillingl (|1977f l. and similarly bv lHavashil (| 198 If ), 
by augmenting the planets' estimated heavy element 
component to solar composition. The necessary mass 
was then distributed in annuli centered on the planets' 
current semimajor axis and a single power law was fit to 
the derived surface density constraints: 
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where £(r) is the surface mass density at a radius r. 

Since its inception, many shortcomings have been rec- 
ognized in the MMSN. Observations of intermediate age 
(2.5 — 30 Myr) clusters indicate a mean disk lifetime of 
~ 6 Myr; consistent with a gas dissip ation timescales for 
circumstellar disks of ~ 1-10 Myr ()Haisch et alj|2001t 
iHernandez et al.ll2007f) . It must be emphasized that the 
MMSN, by definition, contains the minimum amount of 
mass necessary to build the planets at their current semi- 
major axes. Therefore, any proposed solar nebula more 
massive than the MMSN is allowable. Despite this, the 
canonical MMSN has been used for decades as the initial 
conditions for both disk evolution and planet formation 
simulations. It is difficult to grow the cores of the giant 



planets within the time constraint of the gas dissipation 
timescale given the low mass surface densities predicted 
for the MMSN. 

The recent development of the Nice model has shed 
new light on the process of planet formation in the 
solar system (iGomes et all 120051: iMorbidelli et~aU I2005t 
iTsiganis et al.ll2005ft . Themodel assumes the giant plan- 



ets formed in a much more compact configuration than 
they now reside. Simulations suggest that the giant plan- 
ets migrated through scattering interactions with small 
planctesimals exterior to Neptune. The mass and loca- 
tion of this planetesimal disk play critical roles in the 
outcome of their simulations. These interactions caused 
the inward migration of Jupiter and the outward migra- 
tion of the other giant planets. In their simulations, 
chaotic behavior in the outer solar system is initiated 
by the crossing of Jupiter and Saturn through their 2:1 
mean motion resonance (MMR). The chaotic behavior 
causes the rapid outward migration of Uranus and Nep- 
tune. Uranus and Neptune switch places in about half 
of their simulations which nic ely explains why Ne ptune 
is more massive than Uranus ()Tsiganis et alj [2005). 

Assuming the giant planets formed in the compact con- 
figuration predicted by the Nice model, iDeschl (|2007j ) de- 
veloped a new steady state disk model. The predicted 
disk has a much steeper profile and much larger surface 
densities than that predicted by the MMSN. A truncated 
decretion disk, characterized by outward mass transport, 
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is required in order to maintain the steep profile in a 
quasi steady state. Photoevaporation was invoked as a 
natural mechanism for truncating the disk and remov- 
ing mass at the outer boundary. It has been shown that 
such large surface densities can cause significant plane- 
tary migration w ith planets being rapidly lost into the 
Sun (lCridall2009D . 

Star formation in giant molecular clouds implies that 
the majority of stars in our galaxy were born in clusters 
rather than in isolation. The same is most likely true 
for the Sun. ninety percent of stars born in clusters are 
born into rich clusters wi th 100 or more mem bers with 
masses m excess of 50 M^ ([Lada fe Ladall200l . It is esti- 
mated that the Sun formed in a cluster of 1000 — 10, 000 
stars, which in turn implies an average external FUV 
flux that is a few thousand times the non-cluster back- 
ground, but with a stan dard deviation that is compara- 
ble to the average value (Fatuzzo & Adams 2008; Adams 
2010). The compelling reasons for the Sun being formed 
in a cluster of this size are (1) the abundance of the 
short-lived radioactive nuclide 60 Fe derived from mete- 
oric samples cannot be produced by spallation in the 
solar system and can only be explained by an extraso- 
lar nucleosynthetic origin. The capture of this extraso- 
lar 60 Fe into the early solar system is more likely if the 
Sun's birth clust er contained a sufficient number of mas- 
sive supernovae (|Wadhwa et all IIM IAdamsll20Toh . (2) 
Sedna's orbit requires a stellar enc ounter at a distance 
of less than or equal to 400 AU dKenvon fe Bromlevl 
[200llMorbidelli fe Levisodl200l iBrasser et al.ll2006D . A 
smaller birth cluster would not yield such a close stel- 
lar encounter during its lifetime and cannot provide the 
necessary amount of 60 Fe whereas a larger birth cluster 
would give such a large FUV flux that it would pho- 
toevaporate t he disk before the giant planets can form 
(lAdamsll2010l) . 

Observations of low mass young stellar objects near 
the Trapezium cluster in Orion show disks silhouetted 
against the background nebula, the so called proplyds 
(PRoto PLanetary DiskS; Bally et al. 1998, McCullough 
et al. 1995). Initial modeling of these disks invoked pho- 
toevaporation as a result of ionizing radiation from an 
embedded central source. Observations of the same clus- 
ter show that mass loss must be due to neutral flows 
generated at the disks' surfaces. The outflows then be- 
come ionized at some d istance from the base of the flow 
([Johnstone et al.l [19 98) . A natural explanation for the 
observed neutral outflows is that the disks are heated by 
far-ultraviolet (FUV) radiation with energies in the range 
of 6 — 13.6 eV. The neutral outflows then expand and 
are ionized by extreme ultraviolet (EUV) radiation with 
E > 13.6 eV. Models of external irradiation by nearby, 
massive stars were success ful in explaining the o bserva- 
tions of proplyds in Orion ([Johnstone et al.lll998D . 

The mass-loss rate due to photoevaporation by an ex- 
ternal FUV source can be readily derived using a simple, 
spherically symmetric model of outflow from an infinitely 
dense cloud of radius r c . Assuming that gas, of mean 
molecular weight (/i), is driven outward at a constant 
speed, w w , the mass- loss rate is 

M = 4Trr 2 (n)v w n(r), (2) 

where n(r) is the number density of the flow at radius r. 



The column density of the outflow, Ah, is given by 

/>oo 

A H = / n{r)dr. (3) 

Eqn. ([2]) can then be solved for n(r), substituted into 
Eqn. ([3]) and integrated: 

f°° M , M 
A H = / , ,, , dr = — . 4 

The result is then solved for the mass-loss rate, M. The 
mass-loss rate is dependent on the column density of at- 
tenuation, Ah, and proportional to the radius, r. 

M = 47rr c (/x)« w A H . (5) 

The outflow velocity can then be set equal to the sound 
speed in the FUV heated disk "atmosphere" which is 
analogous to the isotherma l atmosphere in the simplified 
model that was derived by lAdams et al.1 <|2004f ) and will 
be presented in ^3] This results in the mass-loss rate of 
Eqn. ([T3|) differing only by the geometric factor J 7 , which 
incidentally is of order unity. 

A visual magnitude of extinction A v , of order unity, 
typically requires the column density A(H) « 5 x 
10 21 cm~ 2 . A column density of roughly 10 21 cm -2 
is a generally accepted value for comple te extinction 
([Johnstone et all 119981; lAdams et all 120041) . If the disk 
atmosphere is hea ted to a temperature of 1000 K 
([Adams et all 120041 ) . then the gravitational radius, be- 
yond which heated gas can escape from the Sun, is about 
100 AU (see Eqn. (11)). Our estimate of the mass-loss 
rate is therefore M w 1.5 X 10~ 7 M Q yr _1 , which is cer- 
tainly enough to clear the solar nebula on the 10 6 yr 
timescale needed to match observations. 

External irradiation has previously been invoked to 
successfully explain the over abundance of noble gases 
in Jupiter. The Galileo spacecraft made in situ measure- 
ments of the composition of Jupiter's atmosphere. It 
measured an enrichment in heavy elements with respect 
to solar values. In particular, it measured the enrich- 
ment of the noble gases Ar, Kr, and Xe, species that 
only condense at temperatures of less that 100 K. In 
light of these findings, it was proposed that Jupiter must 
have formed from low temperature planetesimals in a 
cold environment, well outside of Jupiter's or bit, allowing 
for th e direct condensation of noble gases ([Owen et all 
[19991) . It was then proposed that the noble gases could 
have been delivered by clathrates, allowing for higher 
formation temperatures ([Gautier et all [20011 Deliver- 
ing a large amount of noble gases as clathrates would 
require a large amount of H2O ice. Based on internal 
structure models of Jupiter that constrain the core mass 
to be < 40 M®, it is difficult to explain how sufficient 
amount s of noble gases were d elivered to Jupiter in this 
fashion (IGuillot fe Huesol [20061) . 

Using the e quations of mass loss from lAdams et all 
(pOOl (see p. IGuillot fe Huesol (f2006l ) were able to ex- 
plain the over abundance of noble gases in Jupiter's at- 
mosphere as a result of the loss of hydrogen and helium 
through the photoevaporative process. They used a sim- 
plified one-dimensional evolutionary disk model, in which 
the vertical structure was averaged, to determine the en- 
richment of the solar nebula in heavy elements. The 
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noble gases are assumed to be trapped within solids in 
the cold outer disk while hydrogen and helium are re- 
moved by photoevaporation. As solids migrated inward 
due to gas drag, they are then released as gases in the 
inner disk and incorporated into the giant planets. This 
scenario differs from previous models that required the 
noble gases to be accreted while trapped in solid plan- 
etesimals. Removing the requirement that noble gases be 
delivered in solids further loosens const raints on nebular 
tempe ratures in the outer solar system. iGuillot fc Huesol 
(2006) modeled photoevaporation from both EUV gen- 
erated by an embedded early Sun and by ambient FUV 
generated by neighboring cluster members. They found 
that the scenario involving EUV was problematic due 
to the long timescale involved in the removal of nebular 
gas and the high constant value of EUV flux from an 
"unidentified mechanism" . In the FUV scenario, they 
found that the observed enrichment in Jupiter's atmo- 
sphere could be matched with a wide range of values in 
their free parameters. They also found that for disk at- 
mospheres heated to T onv > 100 K, the disks dissipated 
o n Myr t imesc ales, consistent with observations. 

iDeschl (|2007l ) derived a steady state decretion disk that 
matched surface density constraints derived by assuming 
the four giant planets formed in the compact configura- 
tion of the Nice model. We extend this line of investi- 
gation by modeling a time-dependent disk that suffers 
mass loss by photoevaporation from an external FUV 
source. Radiation from the central star certainly also 
played some role in disk evolution. Recent numerical 
simulations which combine photoevaporation from the 
central star with viscous evolution show that the ma- 
jority of mass-loss is dominated by loss from the outer 
regions of the disk where less ti ghtly bound mater ial can 
easily escape from the system (jGorti et al.l [20091 . The 
mass loss rates from irradiation by the central star are 
generally lower than those due to external irradiation, 
but depending on the assumed FUV flux they can be of 
comparable magnitude. The FUV environment of any 
given star in a young cluster is highly uncertain. Due 
to the incomplete sampling of the stellar initial mass 
function, the variety of cluster sizes and the 1/r 2 de- 
pendence of flux on the radial distance from the center 
of the cluster, where the most massive stars reside, the 
standard deviation of predict ed FUV fields can be com- 
parable to th eir mean values ([Proszkow fc Adams! 120091 : 
Adama l2010l ) . When the highly variable nature of young 
stellar emission is also taken into account, the ratios of 
internal to external FUV flux become even more uncer- 
tain. Given these uncertainties, it is generally not until 
late times that the internal source begins to affect the 
inner disk; o pening an inner gap and rapidly disp ersing 
the system (jClarke et al.1 [200lt iGorti et al.l[2009h . For 
these reasons we have neglected UV irradiation from the 
Sun in this paper. 

It is our goal to explore the claim that a truncat ed disk 
will pr oduce a disk with the steep profile found by IDeschl 
(|2007t) . We use the outputs from our viscously evolving 
disk model to investigate the growth rates for the giant 
planets within an evolving disk. In accordance with the 
Nice model, w e place our g rowing embryos at the same 
orbital radii as IDeschl (|2007l) . Once the growth rates have 
been determined, we then investigate Jupiter's chances 
of survival against migration because it is the most likely 



giant planet to be lost into the Sun. 

Observations of diverse extrasolar planetary systems 
in recent years make it evident that planetary migration 
plays a significant role in the planet formation process. 
Two types of migration have been widely investi gated, 
comm only known as type I and type II migration (jWardl 
Il997f ). Type I migration occurs when spiral density waves 
are launched in a disk from an orbiting body. The density 
waves exert unbalanced inward and outward torques on 
the orbiting body. The sum of these torques is generally 
negative, causing a body undergoing type I migration to 
spiral inward on a relatively short timescale compared 
to the viscous timescale of the disk. Type II migration 
occurs when a body grows sufficiently large to open a 
gap in the disk. The body is then drawn along and its 
orbit decays on the viscous timescale of the disk. Type 
I migration affects smaller bodies and type II migration 
affects larger bodies. The cutoff between the two types 
of migration occurs when the Hill sphere of the orbiting 
body becomes comparable to the scale height of the disk. 
In high-mass disks, runaway migration may occur before 
the planet can open a gap, leading to very rapid orbital 
decay (type III migration). Despite the current uncer- 
tainties in migration theory, it is almost certain that it 
played a role in the processes that formed giant plan- 
ets. Any comprehensive model of planet formation must 
consider the im plications of migration on planet survival 
(|Armitagdl2"oToh . 

This paper is organized in the following manner. In 
$2] we will review Desch's (2007) MMSN model and it's 
implications for planet growth. In fJ3]we will discuss the 
role that photoevaporation plays in the evolution of the 
solar nebula. In £j4]we will present our model of a time- 
dependent protoplanetary disk that is losing mass due to 
photoevaporation by an external massive star. In <j5j fj6] 
and Sj7]we present our results for evolving disks, planet 
growth within those disks and the chances of Jupiter's 
survival against migration. discusses our results and 
lays out plans for future work. Finally, a brief summary 
will be presented in SjH 

2. DESCH MODEL 

The known inadequacies of the M MSN model a nd the 
development of the Nice model led IDeschl (|2007|) to re- 
investigate the primordi al solar nebul a. Within the con- 
text of the Nice model, IDeschl (|2007l ) has developed an 
new MMSN with Jupiter located at 5.45 AU, Saturn at 
8.18 AU, Neptune at 11.5 AU, and Uranus at 14.2 AU. He 
assumes that Uranus and Neptune switched places dur- 
ing the chaotic period following the crossing of Jupiter 
and Saturn through their 2:1 MMR. These four mass 
constraints were combined with mass constraints from 
chondrules, the asteroid belt, and the disk of primordial 
planctesimals laying outside the orbit of Uranus to de- 
velop a single power-law profile for the primordial solar 
nebula: 

, r \ -1 / \ -2.168 

E < r) = 343 (5y (ioau) «™- 2 ' < e > 

where f p is the fraction of the mass of condensible solids 
in planetesimals. 

The surface density profile of Eqn. [6] is not only more 
massive but much steeper than the canonical MMSN 
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(E oc r~ 3 / 2 ). In a steady state thin a-disk, the surface 
density follows the relation E oc 1/v. Using these as- 
sumptions, along with the a- viscosity prescription, would 
imply that S oc T(r) _1 r -3 / 2 . A surface density pro- 
file consistent with the MMSN would therefore imply 
a constant temperature profile throughout the disk. It 
is generally thought that disks are flared due to irra- 
diation from the central star. It was shown that a 
flared disk with an internal radial temperature distri- 

bution of T(r) = 150 r AU ' can produce a spectral en- 
ergy distributi on that is consistent with observations of 
T Tauri stars (jChiang fc Goldreichl [1997T ) . This has led 
many researchers to use a temperature profile of the form 
T(r) oc r~ x / 2 . Again, using the same assumptions about 
the disk, this implies that the surface density profile 
s hould b e £ oc r . 

iDeschl (|2007t ) first investigated whether a viscously 
spreading accretion d isk o f the type studied by 
Lyndcn-Bcll fc Pringld ()1974h would adequately match 
his surface density constraints as well as timing con- 
straints on planet formation. Although the surface 
density constraints could be matched with a viscously 
spreading disk, such a disk evolves too rapidly to satisfy 
the constraints for planet formation. At 10 AU, planetes- 
imals should have formed by 0.03 Myr, but at this time 
he finds densities that are an order of magnitude lower 
than those implied by the augmented mass of Saturn. 
The viscously spreading disk also has trouble matching 
the density profile of the new MMSN (Eqn. ([6])) at small 
radii for ear ly times and at large radii for late times. 
IDeschl (|2007t ) concludes that the various constraints on 
the surface density and timing of the solar nebula are 
best matched with a steady state pro file. 

Following up on his conclusion, IDeschl (|2007l ) re- 
examined the equations for viscously evolving steady 
state disks. For a viscous disk with surface density E 
and viscosity v., the conservation of mass yields 



dE 
dt 



1 dM 



2~Kr dr 

and the conservation of angular momentum yields 
dV 1 rdr \2n dr J 1 



(7) 



(8) 



where Q is the angular frequency. M is the net flow of 
mass through an annulus of the disk with a negative mass 
flux corresponding to inward accretion. Integrating Eqn. 
((U results in 



-M 2 ^ o on 

r CI - r Ei/— = const. 

2tt or 



(9) 



The constant is evaluated by ch oosing an ap propriate 
boundary condition. This is where IDeschl (|2007l ) diverges 
from previous derivations. In the past, the equation has 
been solved by assuming the dominant boundary is the 
inner boundary and the evolution of the disk is governed 
by the mass flux across the inner boundary while the 
outer boundary is allowed to expand indefinitely such 
that angular momentum is conserved. 

Assuming a temperature profile of the form T{r) oc 
r~ 9 , a viscosity of the form v = ac s /tt, and the standard 



a- viscosity prescription which will be discussed in he 
finds a general solution for the surface density profile, 



E(r) 



M 



2>TTv{r ) \r 



"(3/2-9) 



3TTv(r ) \ r 
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(10) 

Instead of using boundary condi tions t o deter mine the 
constants of integration rh and ro,|Desch (2007) matches 
this solution to his derived surface density profile (Eqn. 
([6])). For a steady state disk to be consistent with a sur- 
face density profile E(r) oc r -2168 , the mass flux M < 0. 
This implies that the solar nebula must be a decretion 
disk, with significant mass loss from the outer disk edge, 
rather than an accretion disk. Unlike a standard accre- 
tion disk which only requires sufficient outflow for the 
removal of angular momentum, a decretion disk is char- 
acterized by an significant outward net flow of mass. The 
outward mass flow is driven by photoevaporation which 
requires a constant supply of mass from the inner solar 
system to replenish mass lost at the outer edges. De- 
pending on local conditions, mass loss from the outer 
edge can dominate over accretion onto the central star 
necessitating an outward mass transport throughout the 
majority of the disk. It is also implied that the disk 
must be truncated at an outer edge, r^. Truncation can 
be naturally explained by invoking photoevaporation by 
an external source. 

Although Desch's MMSN is an improvement over the 
original M MSN, it still suffers its own limitations. For 
one, IDeschl (|2007l ) is unable to constrain his disk's outer 
radius to better th an within 30 — 100 AU. A recent pa- 
per by (Crida 2009) also points out that the large surface 
densities present in Desch's model would cause substan- 
tial migration of the giant planets. Using the hydro code 
FARGO in its 2D1D version, their models show Jupiter 
quickly falling into the regime of type III runaway mi- 
gration and rapidly falling into the Sun. Even though 
IDeschl (|2007| ) was able to construct a steady state solu- 
tion, a steady state profile is an oversimplification and 
is inconsistent with a disk eroded by photoevaporation if 
the planetary accretion timescale is a significant fraction 
of the disk lifetime. 

3. PHOTOEVAPORATION 

The presence of early-type stars in the vicinity of the 
solar nebula would have exposed it to FUV radiation. 
FUV radiation would have heated the periphery of the 
disk. Gas heated to sufficient temperatures would then 
have become unbound from the disk. The gravitational 
radius, r g , is defined as the radius at which the sound 
speed of the heated gas, a s , equals the escape speed from 
the system. Gas beyond the gravitational radius will 
escape from the system: 



kT 



(11) 



Here G is the gravitational constant, M* is the mass of 
the central star, k is the Boltzmann constant, and T 
is the temperature of the super-heated atmosphere, or 
what we will refer to as the envelope temperature, T env . 
The gravitational radius is the canonical radius beyond 
which gas heated to a temperature T env will escape from 
the disk. In actuality, gas can escape from the disk at 
radii substantially smaller than r g . 
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In reality, the heated atmosphere has a depth- 
dependent temperature and the heating and resultant 
outflow are complicated processes. Because of these 
complexities, it is useful to employ a simplified model 
with an isothermal atmosphere. Consider a disk irradi- 
ated and heated by external FUV radiation. Depending 
on the strength of the FUV flux, the heated gas will 
reach temperatures in the range 100 K < T < 3000 K 
(| Adams et al.l 120041 ). As the gas heats, it expands gen- 
erating a neutral outflow. The expanding outflow be- 
gins subsonically but becomes supersonic by the time it 
reaches the gravitational radius. This outflow is generally 
isotropic, but the majority of mass loss is dominated by 
mass loss from the outer edge of the disk. The isotropic, 
neutral outflow serves to shield the disk from EUV radi- 
ation that would ionize the disk and heat it to ps 10, 000 
K. A diagram of a disk of radius r d around a star of mass 
M* illuminated by an external sourc e and the subsequent 
resultant outflow was presented by lAdams et al.l (2004) 
and is reproduced in Figure [1] This outflow is, in effect, 
a super-heated atmosphere that can be characterized by 
a single envelope temperature, T cnv . 




Fig. 1. — (Figure taken from I Adams et all <2004fl .1 Schematic 
of a disk with radius around a star of mass M» , illuminated by 
the FUV (and perhaps EUV) radiation from nearby stars of grater 
mass. The disk is inclined so that the top and edge are exposed. 
The disk scale height is at the outer radius r<j. In the subcritical 
regime, where r<j < r g , the bulk of the photoevaporation flow (the 
radial flow) originates from the disk edge, which marks the inner 
boundary. The flow begins subsonically at r^, with speed and 
density n^. The flow accelerates to the sound speed at r s (the sonic 
point), which lies inside the critical escape radius r g . Beyond the 
sonic point, the flow attains a terminal speed and the density falls 
roughly asnocr -2 . Although some material is lost off the top and 
bottom faces of the disk (the vertical flow), its contribution to the 
mass-loss rate is secondary to that from the edges. Nonetheless, 
the polar regions are not evacuated, the star is fully enveloped by 
the circumstellar material, and the incoming FUV radiation will 
be attenuated in all directions. 



Until recently, only mass loss beyond the gravitational 
radius has been considered. Analytic arguments and nu- 
merical experiments have sho wn that gas can be removed 
down to a radius of (0.2 • r g ) (Liffman 2003; Adam s~t al.l 
2004). Although the heated gas at these radii is pre- 
vented from directly escaping from the disk, there exists 
an atmosphere which can extend beyond r g . This at- 
mosphere can be photoevaporated away and a resultant 
outflow will develop. The outflow will behave very much 



like a Parker wind. As mass is lost from the out-flowing 
atmosphere, it is replenished from the disk and mass is 
effectively lost at r < r g . Assuming a temperature of 
2cnv = 1000 K for the heated atmosphere of the solar 
nebula, r g ~ 100 AU. At (0.2 • r g ) = 20 AU, the forma- 
tion of planets would be effected by the photoevaporative 
outflow. For a more thorough d i scussion of subcritical 
mass l oss, see lAdams et al.l (|2004); Hollcnbac h fc Adamsl 
(l200l . 

In an effort to calculate the probability that a solar- 
type star would experience sufficient photoevaporation 
from exter nal irradiation t hat it would affect giant planet 
formation, lAdams et al.l (|2004l) investigated the mass- 
loss rates from circumstellar disks due to external FUV 
radiation. They studied the previously unexplored sub- 
critical regime, where the outer radius of the disk, r d , 
is sma ller than the gravitational radius. Ada ms et al.1 
(2004) used a photodissociation region (PDR) code to 
determine the depth-dependent temperature of the gas 
based on the optical depth, density, and FUV flux. Their 
PDR code also included 46 chemical species and 222 
chemical reactions. The chemistry is critical for deter- 
mining the cooling rate of the gas. For a given radiation 
field strength Go, disk size r d , disk temperature Tj, and 
stellar mass M*, an iterative procedure was used to de- 
termine the density at the base of the flow n d as well as 
the flow speed at the inner boundary. These two quan- 
tities then determine the mass-loss rate. 

In order to u nderstand the result s of their detailed nu- 
merical model. lAdams et al.1 (|2004[ ) developed simple an- 
alytical models for the photoevaporative mass-loss rates 
for cases in which r d , the location of the outer edge of the 
disk, is both inside and outside the gravitational radius. 
These models are characterized by a single temperature 
T cnv . Although the strength of the FUV radiation field 
Go does not directly enter into these equations it speci- 
fies the envelope temperature which determines a unique 
sound speed, a s for the isothermal atmosphere. 

M = G iVc</i>a s r g (^)e X p(-^) r d < r g (12) 

M = 47rJ r (^)o-p T ] JV a s r d r d > r g (13) 

The first equation is for subcritical disks, and the second 
equation is for supercritical disks. Nc is the critical sur- 
face density of the flow and ctfuv is the cross section for 
dust grains interacting with FUV radiation. The dust 
optical depth is given by rpuv — ofuv ' ^ or an 

optical depth of order unity, Cp^y ~ Vjj, where Ah is 
evaluated at the critical density Nc ~ 10~ 21 cm' 2 . The 
factor Go is a constant of order unity used by Adams 
et al. (2004) to match their numerical and analytical 
solutions. It is used in our model to match the mass- 
loss rates for sub- and supercritical disks at a radius of 
r d/r g = 0.25. This is necessary because the mass-loss 
rates are sensitive functions of the strength of the FUV 
radiation field as well as the assumed matching point. 
We have matched the subcritical solution onto the su- 
percritical solution, because the supercritical solution is 
better understood and well constrained. The factor J 7 , 
in the second equation, is the fraction of the solid angle 
subtended by the flow and is ~ 1 because the flow from 
the disk surface and edge merge at a radius between r d 
and 2r d ; creating a nearly spherically symmetric outflow 
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(|Adams et al.ll2004ft . 

We have chosen to model irra diation from an e x terna l 
source because of the success of Uohnstone et all (|1998[ ) 
in modeling proplyds in Ori on as well as t he need for a 
truncated disk discover edbv lDeschl ((20071) . lAd ams et alj 
(2004) have shown that any incident EUV is attenuated 
very rapidly in the disk atmosphere at several disk radii. 
It photoionizes a portion of the disk atmosphere but 
is unable to penetrate deeply and affect the disk itself. 
EUV radiation can perhaps affect disk evolution at late 
stages when it could help to clear the gas on very short 
timescales. Therefore, our research has been focused on 
FUV radiation and as of yet has neglected the effects of 
EUV radiation. 

4. DISK MODEL 

Under the thin disk approximation the continuity 
equation for a small annulus of Ar located at radius r, 
as Ar — > for a disk with surface density £(r, t) and 
viscosity v is 

as d 



-(rT,v r ), 



(14) 



dt dr 

where v r is the net radial velocity of material transported 
via viscous processes in the disk. In the same disk, the 
conservation of angular momentum is described by 

|(Er 2 fi) + i|-(Sr 3 ^ r ) = llfuT^^), (15) 
at r or r or \ ar J 

where fl is the Keplerian angular velocity, — 
(GM/r 3 ) 1 / 2 . 

Since its inception, the a-viscosity prescription has 
been a useful tool for inves tigating the temporal evo - 
lution of circumstellar disks (|Shakura fc Sunvaevlll973| ). 
Using the a-viscosity prescription, the viscosity can be 
parameterized in terms of the local sound speed in the 
disk, c s , and the local Keplerian orbital velocity, £\ C p : 



ac s H 



2 n -l 



(16) 



Self-consistently solving for the viscosity requires solv- 
ing nonlinear differential equations and is beyond the 
scope of the current work. The mechanisms that generate 
viscosity are not very well understood. Using the admit- 
tedly naive alpha prescription requires iteratively solving 
for the vertical and radial temperature structure. A more 
realistic model would also consider radiative transfer and 
cooling, processes regulated by poorly constrained dust 
opacities. 

For simplicity, we assume the viscosity is proportional 
to the radius of the disk such that 



(17) 



where v and Rq are scalings for viscosity and radius. 
This ass umption has been used in the past by many other 
authors (|Clarkdl2007t iHartmann et al.lll998l) . 

The linear dependence of viscosity on radius in our 
model implies that the temperature profile in our disk is 
proportional to r -1 / 2 . In all of our simulations we use a 
temperature profile of the form T(r) = 150-r -1 / 2 K. This 
is cons istent with e arlier works, and in particular with 
that oflDeschl (12007ft . He us es a temperature profile from 
iChiang fc Goldreichl (|1997f ) that is of the form T(r) = 
150 • r~°- 429 K. Evaluating the temperature profile above 



at r = 1 AU, T disk = T(l AU) = 150 K and using the 
alpha prescription for viscosity, we are able to determine 
the viscosity scaling constant. 



<kT, 



disk 



(a») 



(18) 



where fc is the Boltzmann constant. 

At this point, it is useful to define the variables h and 
g, the specific angular momentum and torque. They are 
good variables to use when the viscosity is proportional 
to r as they allow us to transform the viscous disk equa- 
tion into a simpl e, linear differenti al equation with a con- 
stant coefficient (jHartmannl 1 1 998h : 



and 



h = r 2 n = (GMpr) 1 / 2 

i 

g = —2nrY*vr —— = Snvh'E, 
dr 



(19) 



(20) 



where M p is the mass of the primary and G is the grav- 
itational constant. By expressing Eqn. (|15[) in terms of 
h, it can be shown that they obey the relation 



^ = -M 

dh 



(21) 



where M is the outward mass flux in the disk. 

Using the above relation, the continuity equation be- 
comes 

9S ■ 1 dM (22) 



which reduces to 



dt 2nr dh 

dg = 3v(GM p ) 2 d 2 g 
dt 4h 2 dh 2 ' 



(23) 



By substituting in the functional form of viscosity, 
Eqn. (|17[) . the viscous disk equation can be written as 



dg _ 3 u GM p d 2 g 
~di ~ I Ra dh 2 ' 



(24) 



We are further able to non-dimensionalize the problem 
resulting in a linear first order differential equation, 



where 



dg_ = (Pg_ 

dt ~ dh 2 ' 



(GMptfo) 1 / 2 ' 



and 



t = 



4z/ t 
3R 2 



(25) 



(26) 



(27) 



(28) 



are the non-dimensional variables. 

As with most previous disk models, we employ a zero 
torque inner boundary condition. This allows for accre- 
tion from the disk onto the Sun. The outflowing material 
carries some finite amount of angular momentum away 
from the system and must therefore exert a torque on 
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the outer edge of the disk. The torque exerted on the 
outer edge of the disk can be expressed as 



,9d = 3\/27riV c (^)^o( -^r )<m 



(29) 



We use this torque as the outer boundary condition in our 
models. In terms of our non-dimensionalized variables, 
the outer boundary condition is 



9d 



hi 



(30) 



It must be noted that these two boundary conditions 
are used to solve the equation governing the temporal 
evolution of the mass surface density. We have, in ef- 
fect, a second outer boundary condition that governs the 
temporal evolution of the location of the outer boundary. 
This condition is set by a differential equation that was 
derived using mass conservation at the outer boundary. 



Mbo 



M 



M 



(31) 



The mass flux due to viscous processes is simply 
~dg/dh and the mass flux due to photoevaporation is 
taken to be Eqn. (fT2")) or Eqn. (fl"3")) depending on the 
location of r&. The mass flux due to the motion of the 
boundary is 2itrd^d(drd/ dt). The equations governing 
the location of the outer boundary was derived by sub- 
stituting these expressions into Eqn. (f3"l"|) . Depending 
on the location of the outer boundary, we will either be 
in the subcritical regime, 



1 



dg C a s r„ 



dd = 

'H V2nN c ((i) dh d 

or the supercritical regime, 





-Is 


[ — ) exp 









dt 



1 



dg 



2ttN c {^jl) dh d 



- v2~7r a s r d . 



(32) 



(33) 



Expressed in terms of our non-dimensional variables, 
the motion of the outer boundary in the subcritical 
regime is 

dhd 
dt 



1 



hi 



• s 



1/2 



at. 

dhd 



C a s r„ 



V h\ 

and in the supercritical regime is 



exp 



difference scheme, 
At 



and 



im+l — 



dhd 
dt 



1 



1/2 



dg 
dhd 



87r a s i?o 
3 v a hd 



2R h 2 d 
(34) 



(35) 



Using th e formalism of t he va riable space grid (VSG) 
method of iKutluav et al.l (|1997f ) (see the Appendix for 
details), Eqn. (p5|) can be discretized with an explicit 
finite difference method: 



2 Ah si \R 



1/2 



-4^_!+^_ 2 )-At 



~2 C a s r„ 



3v8tt vq Rq 
(37) 



At 



2 Ah si 



1/2 



Tt m-^_ 1+ ~g^_ 2 )-At 



8n a s R a 
3f s m 



(38) 

At each time step, depending on whether we are in the 
sub- or supercritical regime, either Eqn. (|3"T)l or Eqn. 
(|38p must first be solved to determine the new boundary 
location. Then the viscous evolution, Eqn. is solved. 
All of the simulations presented here were performed on 
a Mac G5 PowerPC. They were evolved on a grid of 200 
points evenly spaced in specific angular momentum. 

Our i nitial conditions are that of t he similarity solu- 
tions of iLvnden-Bell fc Pringli (| 19741 ). 



Mo 
2irR 1 r 



exp 



r 

R~i 



(39) 



where R\ is the initial disk scaling radius and Mq is the 
initial disk mass. 

5. DISK EVOLUTION 

The two free parameters of our model that most di- 
rectly influence the timescales of disk evolution and 
planet formation are the temperature of the super- heated 
atmosphere, T onv , and the non-dimensional viscosity pa- 
rameter a. These are also the two parameters that are 
the least constrained. We have performed a number of 
simulations to investigate the role played by these two 
parameters in affecting disk morphology and evolution 
timescales. The various input parameters and resul- 
tant timescales are tabulated in Table Q] Our reference 
model best matches the parameters, temperature, vis- 
cosity, etc., of Desch's model and canonical disk models 
in general. The reference model has an envelope tem- 
perature of T env = 600 K and a viscosity parameter of 
a = 0.001. 



TABLE 1 

Input parameters for various models. 



simulation 


T cnv [Kj 


a 


normalized timescale 


LV 


600 


0.0001 


3.6 


LT 


100 


0.001 


3.6 


Reference 


600 


0.001 


1.0 


HT 


3000 


0.001 


0.45 


HV 


600 


0.01 


0.81 



- r+ l = - r . 



At 



2 s m Ah 



+i 



(36) 

where r is defined as At/(Ah) 2 . The equations govern- 
ing the location of the boundary (Eqn.'s (f34|) and (|35|) ) 
can then also be discretized with the same explicit finite 



In addition to our reference model, a number of sim- 
ulations have been completed in order to investigate the 
effect that the adjustment of key parameters can have 
on the evolution of the disk and on planet formation 
9™- 1 ) + r (<?f+ 1 - 2 dT +<?!- i ) timescales within that disk. For these simulations, we 

have varied the envelope temperature and viscosity pa- 
rameter a to their extreme values as predicted by canon- 
ical disk and complex PDR models. 

The viscosity parameter a in protoplanetary disks 
is considered to lie within an order of magnitude of 
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0.001. We have therefore considered values of a be- 
tween 0.01 and 0.0001. According to complex PDR 
models that have been used to determine temperature 
and density profiles of photoevaporating outflows, the 
envelope temperatu re lies in the range 100 — 3000 K 
(|Ad ams et~ayl200l . The given range of temperatures 
corresponds t o FUV fields with 300 < G < 3000 
(I Adams et alJ feOO^). G is a dimensionless quantity ex- 
pressed in terms of the Habing field (1 Habing field = 
1.2 x 10~ 4 erg cm" 2 s" 1 sr _1 ) and where Go = 
1.7 Habing fi elds for the local interstellar FUV field 
(|Tielensll2005D . In an effort to constrain the behavior of 
our photoevaporative disks, we have used only extremes 
in these values. The models have been labeled according 
to which parameter (s) have been varied, with L or H re- 
ferring to either low or high values of either the viscosity 
(V) or the temperature (T). 

From these simulations, we have determined the time 
evolution of the surface density at the locations of the 
four giant planets in the compact configuration of the 
Nice model. Using the time-dependent surface densities, 
we then estimate the growth rates of the giant planets' 
cores. The growth rates and decaying surface density 
profiles were then used to calculate the migration rates 
of the giant planets. 

In general, the evolution of the disks produced by 
our numerical simulations begin with a prolonged phase 
where the outer boundary expands as viscosity trans- 
ports mass outward from the massive inner disk. As the 
disk evolution continues, it is characterized by the slow 
erosion of the outer boundary and a nearly self similar 
shape of the disk's radial mass surface density profile. 
One feature all disks have in common is a mass front 
located at a truncated outer boundary. 

With the exception of the reference model, the models 
were run until the mass, M comp , of the gas disk within 
the region of giant planet formation, 2 AU < r < 30 AU, 
had reached a given value. We have chosen this region 
because it al lows us to compare our models to that of 
iDeschl (|2007l) . His steady state disk model contains ~ 
0.07 M© in the comparison region. Our models were 
run until M comp = 0.07 M , 0.035 M Q and 0.0175 M©. 
The reference model was additionally run until Af comp = 
0.00875 M . It took the our reference model 0.70 Myr 
to evolve from a mass of 0.07 M© to a mass of 0.035 M© 
and 0.38 Myr to evolve from there to a disk mass of 
0.0175 M©. Throughout this work, all times listed for a 
given specified model are relative to the time when that 
particular model contains M comp = 0.07 M©. 

Snapshots of the evolving mass surface density of the 
reference model are shown in Figure [2] at four different 
times that span a 1. 3 Myr i nterva l. The surface den- 
sity constraints from IDeschl (|2007f) . as well as his de- 
rived surface density profile, have been overplotted for 
comparison. The uncertainty in the solid component of 
the two inner giant planets, with Jupiter in particular, 
is large because of the inability of hydrostatic models, 
which rely on high-pressure equations of state, to con- 
strain current core masses. By inspection, one can see 
that as the disk e volves i t mat ches with the surface den- 
sity constraints of[fc)csch (2003) at various radii at various 
times. The inner surface density constraints are matched 
early on and the outer constraints are matched at later 



times. The ad hoc profile of Eqn. is plotted with 
the dotted line. It is interesting to note that although 
our first output profile contains the same amount of mass 
as his profile between 2 and 30 AU, his profile is much 
steeper than our model profiles. Our surface density pro- 
files have £(r) oc r l zo -o.33 5 where the range of expo- 
nent here is not the uncertainty but represents rather 
a range of profile slopes. All profiles were fit with a 
power-law slope through the giant planet-forming region, 
5 AU < r < 15 AU. 

Generally speaking, we have two families of disks, those 
with contracted radial surface density profiles and those 
with extended profiles. The contracted disks are domi- 
nated by photoevaporation and have slopes with an aver- 
age power-law slope of —1.6, whereas the extended pro- 
files are dominated by viscous spreading and have an av- 
erage slope of —0.94. The models with contracted pro- 
files and those with extended profiles both have slopes 
that slightly increase with time as the outer disk radius 
moves inward. The steepest profile, and our most con- 
tracted (rd ~ 19 AU), is from the simulation LV and 
has a radial surface densi ty profi l e of r ~ 21 , which is as 
steep as that derived by IDeschl (|2007l ). Although the 
slope of the radial surface density profil e of simulation 
LV matches that derived by Desch (2007), it matches at 
only very late times and does not maintain a steep slope 
such as implied by Desch's quasi steady state model. The 
differences between most o f the d e rived surface density 
profiles of this work and of IDeschl ([2007) probably arise 
from the different assumptions that were made in these 
models. This will be discussed further in Sj8] 




radius [AU] 

Fig. 2. — The radial mass surface density at four times for the 
reference model. Outputs for all simulations look very similar these 
profiles due to our constraint on the mass contained within the 
planet forming region. These surface densities evolved over a 1.3 
Myr time interval. The times plotted, relative to our first output, 
are as follows; 0.70 Myr, 1.1 Myr, and 1.3 Myr. The surface density 
const raints inferred from the Nice model are over plotted JDesch] 
120071) . The dotted line is the surface density derived by IDeschl 
120071) . 



The initial mass is different for each of our models. 
This was necessary such that each model exhibited the 
same behavior at the times of interest. This has however 
required us to use some disks with unrealistically large 
masses, some as large as a 0.5 M©. Such large disks 
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would likely be susceptible to gravitational instabilities. 
We have performed an additional simulation with a small 
initial disk mass (Mo =0.1 M©) so that the behavior 
exhibited by the high-mass models can be verified under 
more physically realistic conditions. A plot of the radial 
surface density profiles at various times (times shown in 
inset) is shown in Figure [3J The elapsed time between 
the first and last outputs is 2.0 Myr. As with all of our 
other models, this mod e l uses the similarity solutions of 
iLvnden-Bell &: Pringld (| 19741 ) as the initial conditions. 
In this case, the initial disk mass is 0.1 M©and the scaling 
radius, Ri, has been set to 10 AU. 




100 
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radius [AU] 

Fig. 3. — The radial mass surface density at five times for the a 
model with an initial disk mass of 0.1 Mq. These surface densities 
evolved over a 2.0 Myr time interval. The times plotted, relative 
to our first output, are show n in the inset box. This mod el also 
uses the similarity solutions of Lyndcn-Bcll & Pringlc ( 1974) as the 
initial conditions. 



Again, as with our other models, the disk begins 
with a rapid expansion of the outer boundary. As the 
disk spreads the slope of its radial surface density pro- 
file quickly approaches a power law with an exponent 
—1.05. This is similar to the average slope from all 
other models presented here. At t ~ 6.6 x 10 5 yr the 
outer boundary of the disk reaches its maximum value 
then begins the shrink. The disk then shrinks, both in ra- 
dius and in overall magnitude, as the outer edge is eaten 
away by photoevaporation. The disk maintains a nearly 
self-similar shape until the outer boundary shrinks con- 
siderably at which point the slope of the radial surface 
density profile steepens slightly. At the end of this sim- 
ulation the disk contains 0.0035 M©. We infer from this 
model that the behavior seen in our models is indepen- 
dent of disk mass. 

Simulations LT and HT were designed to test the effect 
of FUV radiation on the evolution of the disk and in turn 
on planetary formation timescales. There are few, if any, 
strong constraints on the FUV environment of the early 
solar system. We therefore explored a range of disk enve- 
lope temperatures as the envelope temperature certainly 
affects the timescale for disk evolution. It was expected 
that a larger FUV flux and hence a larger envelope tem- 
perature would cause the disk to evolve more quickly. 
Furthermore, it was believed that a lower FUV flux and 
envelope temperature would result in a prolonged disk 



evolution allowing more time for the outermost giant 
planet cores to form. 

Compared to the reference model, the evolutionary 
timescale was in fact larger in LT and smaller in HT. 
It took the disk in LT 3.6 times as long to evolve from a 
mass of 0.07 M© to a mass of 0.0175 M© as the reference 
model disk, whereas the disk evolution in HT was faster 
than that in the reference model by a factor of 0.45. 

Given the wide range of a generally used in solar neb- 
ula models, we varied the value of a to see how it affects 
planetary growth timescales. The viscosity parameter a 
was varied to a value of 0.0001 for LV and to a value 
of 0.01 for HV. As expected, the simulations evolved 
more rapidly for increasing values of a. The evolution- 
ary timescale of LV was 3.6 times longer than in our 
reference model and the timescale of HV was 0.81 times 
shorter than that of the reference model. The viscosity 
parameter has been changed by an order of magnitude 
and shows that, at the reference temperature, the disk's 
temporal evolution is just as dependent on the viscosity 
as the temperature. 

It is interesting to note that while LV and LT both pro- 
duce longer evolutionary timescales than the reference 
model, the surface density profiles generated by them 
are strikingly different. The low viscosity of model LV 
prevents the mass in the inner regions of the disk from 
spreading outward and maintains the outer edge of the 
disk within a few times 10 AU. In contrast, the rela- 
tively high reference viscosity (a = 0.001) of LT allows 
for the massive inner disk to rapidly spread outward to 
> 100 AU where it is slowly eroded by the photoevapora- 
tive outflow. The radial profiles of LV and LT can be seen 
in Figure [4J along with the reference model for compari- 
son. Each disk radial surface density profile corresponds 
to a time when the mass contained within the planet 
forming region, 2 AU < r < 30 AU, is 0.035 M©. These 
correspond to times of 2.0 Myr, 2.1 Myr, and 0.70 Myr for 
LV, LT, and the reference model respectively. Here one 
can clearly see the difference in the radial distribution 
between these two models. A similar dichotomy is seen 
with simulations HV and HT. They both produce short 
evolutionary timescales relative to our reference model, 
but are opposite with regards to the radial extent of the 
surface density profiles they produce. 

6. EMBRYO GROWTH 

If one assumes that planets form via the accumula- 
tion of smaller bodies and not through direct gravita- 
tional collapse, the early stage of planetesimal accre- 
tion is characterized by a period of runaway growth 
(|Wetherill fc Stewartfe 'SQI. During runaway growth the 
velocity distribution of planetesimals is dominated by in- 
teractions with other planetesimals. During this time the 
velocity dispersion of planetesimals is low and gravita- 
tional focusing is effective. While gravitational focusing 
is effective, the largest bodies grow much more rapidly 
than smaller bodies and a bimodal size distribution is 
achieved. 

This is followed by a phase of oligarchic growth, where 
the velocity distributions are dominated by interactions 
of the larger planetary embryos. During oligarchic 
growth the presence of large bodies enhances the velocity 
dispersions of smaller bodies and decreases the velocity 
dispersion of the largest bodies. The increased disper- 
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Fig. 4. — Surface density profiles of the extended disk for LT, 
LV and the reference model. The reference model is shown in 
solid, LT with a dotted line and LV with a dashed line. Each disk 
radial surface density profile corresponds to a time when the mass 
contained within the planet forming region, 2 AU < r < 30 AU, 
is 0.035 Mq. These correspond to times of 2.0 Myr, 2.1 Myr, and 
0.70 Myr for LV, LT, and the reference model respectively. 



sion in velocities of the smaller planetesimals decreases 
the effect of gravitational focusing and the largest bodies 
begin to decrease their growth rate. The system becomes 
dominated by a few large bodies, an oligarchy, separated 
by a few mutual Hill radii. Since our simulations take 
place while large amounts of gas are present, we only 
consider runaway growth. 

Early analytical models, by iSafronovl ([1969ft and oth- 
ers, overestimated the growth timescale of planets by 
upwards of 5 orders of magnitude and were inconsis- 
tent with observational constraints of protoplanetary sys- 
tems tha t show ed the removal of gas in ~ 5 — 10 Myr. 
iLissauerl (|1987l ) developed an analytic model for the run- 
away g rowth of planetary embryos (Eqn. 3 from ILissauerl 
(1987)). To evaluate the growth rates of th e giant planet 
embryos, we use Eqn. (14) of ILissauerl (Il99l . The 
growth rate is defined for an embryo of mass M e 
and radius R e and escape velocity i> osc embedded in a 
swarm of planetesimals with a local surface density S p 
and velocity dispersion a: 



dM c 
dt 



V3 



CSC 



(40) 



The numerical prefactor depends on the velocity distribu- 
tion of planetesimals and many values have been quoted 
in the literature, the value used here of v / 3/2 is due to 
an isotropic velocity distribution. We make the conserva- 
tive assumption that the surface density of solids, E p , is 
0.014 times the surface mass density of gas and that the 
solids-to-gas ratio does not change with time or radius. 

One can see from Eqn. [40] that the growth rate is 
dependent on the geometric radius of the embryo, ttR^, 

2 

enhanced by a gravitational focusing factor, (1 + -*¥■). 
The exact value of the gravitational focusing factor has 
been the subject of much study over the years and is still 
much debated. It has been studied with both analytical 
and numerical studies in a variety of different regimes 
including gas-free and gas-damped accretion. 



Numerical experiments show that the eccentricity and 
inclination of the planetesimals in a swa rm are damped 
due to the i nteractions with gas in a disk (Ko kubo fc Ida! 
[I<M Il998l I2CJ001 12001 . The damping of inclination and 
eccentricity due to gas drag causes, at least at the small 
end of the size distribution, the planetesimals to be in the 
shear-domi nated regime where gravitational focusing is 
important ( Rafikovll200l . We can define a characteristic 
velocity 



I'H = 



GM C 



based on the definition of the Hill radius, 



3Mr. 



1/3 



(41) 



(42) 



where a is the semi-major axis of the embryo. Our char- 
acteristic velocity marks a transition between the shear- 
dominated and dispersion-dominated regimes. When the 
velocity dispersion of planetesimals is smaller than our 
characteristic velocity, a < v u accretion wil l proceed in 
the shear-dominated regime (|Armitagdl2010t l . We adopt 
the characteristic velocity, i>h, for the value of the veloc- 
ity dispersion of planetesimals, er, in all of our calcula- 
tions of embryo growth. When a < Uh, the system is 
in the shear dominated regime and three-body dynamics 
become important. Therefore, Eqn. (|40|) is not strictly 
valid as it is derived considering only two-body effects. 
Also, as the embryos grow the system will transition to 
a dispersion-dominated regime where the embryos will 
grow in an oligarchic fashion. Due to the uncertainties 
in when this transition occurs we have focused only on 
runaway growth. Owing to the large variation in esti- 
mates planetary formation timescales and the wide array 
of unknown parameters, disk mass, viscosity, gas/solid 
ratio, etc., our calculations are not meant to definitively 
describe embryo growth but to illustrate how various sur- 
face density profiles determined by viscous evolution and 
photoevaporation effect planet growth. In this regard, 
the following results on embryo growth should be treated 
with some caution. 

In order to determine the cores' growth rates, it was 
necessary to determine the time evolution of the mass 
surface density. The temporal evolution of surface den- 
sities at the location of each core were fit with decaying 
exponentials. This seemed to give a good fit to the data. 
These fits were then used in Eqn. (|40"| to determine the 
masses of the giant planet embryos as functions of time. 
It should be noted that in our models it is the time- 
dependent surface density that determines the growth 
rates and hence the embryo masses. In our models, the 
embryos initially grow rapidly because the surface den- 
sities are large, but the growth rate then begins to wane 
because the disk evolves and the surface densities become 
small. This differs from most models in which the growth 
rates are large throughout the duration of the simulations 
because they assume unrealistically large, steady state 
surface densities. Our research indicates that, because 
of the similarity of relevant timescales, planet formation 
models must take into account the time-dependent be- 
havior of the solar nebula. 

Embryo masses for the four giant planets as a functions 
of time, determined for our reference model, are shown 
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in Figure O We were able to quite easily build the core of 
Jupiter to > 10 M^well within the ~ 5 — 10 Myr window 
implied by gas dissipation. One can see from Figure [5] 
that, in our reference model, we are unable to grow the 
cores of Saturn, Neptune and Uranus to 10 M e within the 
^5—10 Myr time constraint. The gas simply diss ipated 
too qu ickly for them to form in the allotted time. iDeschl 
( 2007) was able to grow cores of sufficient mass in his 
models because he relied on steady state models with 
surface mass densities that remained large throughout 
the planetary growth process. We feel that our decaying, 
time-dependent model is a more realistic representation 
of the solar nebula. 




2 4 6 8 

time [Myr] 



Fig. 5. — Growth of giant planet embryos in our reference model. 
The masses of planet embryos, plotted from top to bottom, are 
Jupiter, Saturn, Neptune and Uranus. 



Figure |6] shows the growth of planetary embryos in sim- 
ulation LT; a disk with a heated envelope temperature of 
100 K. Because of a smaller mass loss rate at the outer 
boundary, the evolutionary timescale of LT is a factor 
of 3.6 over the evolutionary timescale for the disk in the 
reference model. In this model the cores of Jupiter and 
Saturn were both able to grow cores of 10 M^or more 
during the first 10 Myr of evolution. This is not surpris- 
ing considering the prolonged temporal evolution of the 
disk in the low-temperature model. Despite the success 
at the growth of the cores of two innermost giant plan- 
ets, the cores of Neptune and Uranus are unable to grow 
large enough during the duration of the simulation. 

Figure [7] shows the growth of planetary embryos for 
LV (a — 0.0001). The embryos in this model grow faster 
than the embryos in the reference model, and similarly 
to LT. All embryos, with the exception of that of the 
outermost giant planet, are able to grow to masses of 10 
Mq within the allocated ^5 — 10 Myr. As seen before in 
the low-temperature model, the prolonged temporal evo- 
lution of the low-viscosity model provided a sufficiently 
high surface mass density for a long enough time for the 
three innermost cores to grow to sufficient masses. The 
effect of varying the viscosity has nearly the same effect 
on the evolutionary timescale as in the above case where 
T cnv was varied, but the truncated disk in LV provides 
more mass in the giant planet forming region. The extra 
mass provides a more conducive environment for embryo 
formation as seen with the success in growing the core of 
Neptune to > 10 M e . 

As stated earlier, the surface density of solids, E p , was 
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time [Myr] 

Fig. 6. — Growth of giant planet embryos in LT, where the 
envelope temperature, T en v, has been reduced to 100 K. The masses 
of planet embryos, plotted from top to bottom, are Jupiter, Saturn, 
Neptune and Uranus. 
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Fig. 7. — Growth of giant planet embryos in LV, where the 
viscosity parameter, a, has been reduced 0.0001. The masses of 
planet embryos, plotted from top to bottom, are Jupiter, Saturn, 
Neptune and Uranus. 



assumed to be 0.014 times the surface mass density of 
gas. This estimate is based on the canonical gas/solid 
ratio of 70 derived fro m composition of Comet Halley 
(|Jessberger et alJ ll988). It should be noted that this 
estimate is based solely on the content of H2O ice. Ob- 
servations of the ejecta of Comet 9P/Tempcl 1 during 
Deep Impact showed signifi cant amounts of CO, CO2 
and CH3OH (jA'Hearnl 12008). These ices would certainly 
be present at the locations of Neptune and Uranus and 
would result in a higher solid/gas ratio. Combined mod- 
els of viscous disk evolution and kinetic ice formation 
show an increase in the solid/gas ratio with radius. By 
following a chemical reaction network tracing the forma- 
tion and freeze out of ices in a viscously evolving disk, it 
was found that the solid surface density at Saturn (9.5 
AU) is roughly three times that used in previous mod- 
els of planet formation and that the solid surface den- 
sity at Uranus (20 AU) was highe r by a factor of nearly 
4.5 (jDodson- Robinson et al.l[2009t ). An increase in solid 
surface density can also facilitate the formation of plan- 
ets in a more dramatic fashion. The increase in solids 
would certainly decrease the formation timescale of the 
outermost giant planets. Settling of dust to the disk mid- 
plane and preferential photoevaporation of gas can lead 
to a significant increase in the dust-to-gas surface density 
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ratio. This increase in solid surface density can poten- 
tially become unstable to gravita tional collapse and trig - 
ger rapid planetesimal formation ()Throop fc B ally 2005). 

7. MIGRATION TIMESCALES 

Viscous torques from density waves launched in a 
disk fro m an orbitin g plan et are thought to cause mi- 
gration (iWardl 119971 Il998f) . Using numerical results, 
iTanakaet all (|2002Tl were able to constrain analytical 
models for the torque exerted by corotational and Lind- 
blad resonances on a body orbiting in an isothermal disk. 
The net torque in three-dimensions is 



T = (1.364+0.541— -) 

aa c 



M e a c Q c 



Zeajnl (43) 



where the subscript "e" indicates the values of these vari- 
able at location of the embryo. Here, <z e refers to the em- 
bryo's semimajor axis and c s is the local sound speed of 
the disk. The local orbital velocity fi e is approximated 
by the Keplerian orbital velocity fikcp- Eqn. (|43p can 
then be used to determine the type I migration timescale 
using 

Le = M c {GM Q R e y/ 2 
2T 2T 



(44) 



Figure |8] shows Jupiter's migration timescales for LT 
and HT, models where the envelope temperature, T cnv , 
has been varied along with our reference model for com- 
parison. They have been calculated keeping the semi- 
major axis of Jupiter fixed at 5.45 AU. In general, the 
simulations that grew planets the fastest also suffered 
the shortest migration timescales. At early times, the 
migration timescale decreases due to the growth of the 
planet, but at late times the decaying surface density of 
the disk causes the migration timescale to increase. In 
our reference model, it is unlikely that Jupiter would sur- 
vive orbital decay due to type I migration. However, if it 
can survive through the period in which type I migration 
timescales reach a minimum it has a chance for long term 
survival. 




time [Myr] 

Fig. 8. — Jupiter's migration timescales for our reference model 
and models where the envelope temperature has been varied. From 
top to bottom, the migration rates shown are from HT, the refer- 
ence model and LT. 



Figure [S] shows that the migration timescales of Jupiter 
for the reference model are intermediate to those in LT 
and HT. The migration rates in these simulations are 



mainly affected by the large timescale variations in the 
evolution of the surface density of the disk. The long 
(short) timescales produced by lowering (raising) the 
disk envelope temperature allow planets to grow faster 
(slower) and maintain surface densities at higher (lower) 
levels. The combined effect of larger (smaller) planets 
and higher (lower) surface densities combines to cause 
shorter (longer) migration timescales than in our ref- 
erence model. The migration rates of LV and HV are 
similar to those in LT and HT. Figure |9] again shows 
that models with short evolutionary timescales produce 
smaller embryos in a less massive disk and therefore these 
planetary cores have longer migration timescales. 
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Fig. 9. — Jupiter's migration timescales for our reference model 
and models where the viscosity parameter a has been varied. From 
top to bottom, the migration rates shown are from HV, the refer- 
ence model and LV. 



The migration rates derived using Eqn. (|44l) are for 
type I migration and are only valid when the Hill radius 
of an embryo is smaller than the scale height of the disk. 
When the planet mass exceeds some critical value, the 
migration switches to type II migration and the evolution 
of the embryo's semimajor axis becomes locked into the 
viscous evolution of the disk. Our disk, with its large 
amount of outward transport of material and truncated 
outer radius, could cause an embryo to migrate cither 
inward or out ward depend ing on the semimajor axis of a 
given embryo (Ward 2003) . 

In a truncated disk with outward mass transport, there 
is a critical radius inward of which the mass moves inward 
and is accreted onto the central object and outward of 
which the material is transported outward and eventually 
out of the system. The survival of a growing embryo 
against the effects of type II migration depends on which 
side of the critical radius it is. To investigate the location 
of the critical radius in our model we have calculated the 
radial velocity of the material in our disk: 



i) 



E rV2 dr 



[vH r 



1/21 



which in terms of our non-dimensional variables is 



f'r 



3^o 
Ro 



h_dg_ 

2gdh 



(45) 



(46) 



The results of these calculations for the reference model 
can be seen in Figure [TO] We have plotted the radial 
velocity of the flow for the four times shown in Figure [2j 
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The critical radius lies at 38 AU, 25 AU, 18 AU, and 14 
AU at the four times shown. The times shown occur at 
Myr, 0.70 Myr, 1.1 Myr, and 1.3 Myr, respectively. The 
critical radius, where the curves intersect v r = 0, moves 
inward with time. The critical radius at the times shown 
is exterior to the orbit of Jupiter, but it is rapidly mov- 
ing inward and will at some point transition to a radius 
smaller than Jupiter's semimajor axis. If this happens 
early enough it could save Jupiter from migrating into 
the Sun. This could also affect the radial diffusion of 
dust pa sk. 
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Fig. 10. — The radial velocity of material in the disk in the 
reference model. The four times plotted in Figure[2]have also been 
plotted. They are sequential in time from right to left. The critical 
radius lies at 38 AU, 25 AU, 18 AU, and 14 AU at the four times 
shown. The times shown occur at Myr, 0.70 Myr, 1.1 Myr and 
1.3 Myr respectively. 



8. DISCUSSION 

The five simulations presented here were designed to 
have comparable masses throughout the planet-forming 
region (2 AU < r < 30 AU). Although their radial sur- 
face density profiles are nearly identical within the giant 
planet-forming region, their profiles are very different in 
the outer regions of the disk. The various disk models 
have very different radial profiles at large radii, but one 
must keep in mind that there is in actuality very little 
mass in these regions. Despite the very low mass sur- 
face densities at large radii, these disks can extend, in 
some cases, many hundred of AU from the Sun and may 
have strong implications for the development and evo- 
lution of the Kuipe r Belt and trans-Neptunian objects 
(| Adams et al.ll2004| ). 

Despite many of our simulations having large outer 
disk radii, two simulations, HT and LV have outer disk 
radii of roughly 20 — 30 AU. This is within the current 
semi-major axis of Neptune. The simulation HT has a 
disk evolution timescale that is much shorter than our 
reference model and fails to grow the giant planets within 
the timescale of gas dissipation. In contrast, LV had a 
longer disk evolution timescale than our reference model 
and was more effective at growing the giant planet em- 
bryos. The combination of low viscosity and relatively 
high FUV flux prevented the gas from spreading out- 
ward and created a steeper more compact disk. These 
results show that photoevaporation can affect the solar 
nebula in the giant planet formation region. 

We have run our models assuming a constant external 
FUV flux. There are compelling reasons to believe that 
the incident flux is not constant. Individual stars in the 



solar birth cluster would certainly have experienced mo- 
tion relative to one another. Relative motion between the 
Sun and any of its high-mass bret hren would certainly 
have caused the F UV flux to vary (jProszkow &; Adaing 
amsll2010| ). The observed spread in stellar ages 
within young clusters is usually ~ 1 Myr, but subgroups 
within the same cluster have been observed to have a 
roughly 10 Myr difference in ages d Jeffries et alJl2006h . 
This would imply a highly varying UV environment as 
new stars are born into clusters with the early-type stars 
rapidly moving onto the main sequence on timescales of 
10 4 - 10 5 yr. 

The lifetimes of early-type stars are comparable to the 
lifetimes of circumstellar disks and therefore any changes 
in luminosity experienced during their short lives will 
affect the local UV environment. The FUV flux from 
early-type stars can vary substantially on very short 
timescales (10 4 yr) as they transition through the lu- 
minous blue variable (LBV) stage. LBV stars undergo a 
phase marked by high mass loss and instability. Dur- 
ing the LBV phase, a constant bolometric luminosity 
(L « 10 5 — 10 6 Lq) is maintained, but the dense stel- 
lar wind absorbs EUV such that the majority of the flux 
escapes in the FUV. Observations of proplyds in the Ca- 
rina Nebula (NGC 3372) suggest these outbursts can 
have dramatic conseque nces for nearby protoplanetary 
disks ([Smith et al.l 120031 ). This would imply direct con- 
sequences on disk evolution and survival. In the future, 
it would be interesting to model a variety of FUV irradi- 
ation scenarios and investigate the effect of variable flux 
rates. 

In our models, a was held constant throughout the en- 
tire disk and v oc r which follows from T cx r" 1 / 2 . In 
reality, v should not be a simple monotonic function if 
the ionization fraction of the disk varies with radius lead- 
ing to dead zones where the magnetorotational instability 
is suppressed. We plan, in the future, to model viscos- 
ity in a more self-consistent manner by utilizing a disk 
model with an evolving radial temperature profile. This 
would at first involve using the a-viscosity prescription 
with constant a, but with viscosity v that is a function of 
the local disk temperature. By allowing the viscosity to 
vary throughout the disk it could substanti ally alter the 
rates of planet formati on in various regions dKretke et al.l 
2009; Ly ra et al.1 12009). Without further investigation, it 
is hard to say exactly what the effect of allowing the vis- 
cosity to depend on the local qualities of the disk would 
have on the growth rates and chances for survival against 
migration of growing embryos. 

Furthermore, our use of a constant a could explain the 
discrepancy found between the slope of the radial surface 
densit y profile of our models and that derived by IDeschl 
(|200l . If a was allowed to be lower in the inner portion 
of the disk than in the outer regions, then the shallow 
profile that our models have produced would probably 
not develop. In t his case, it is possible that the steep 
profile derived by IDeschl (|2007| ) would develop. At this 
time, our numerical model is unable to handle a change 
of a in the inner portion of the disk and accommodating 
such a change would require significant modification to 
our code. We think this is an interesting proposition and 
would like to investigate it further in the future. Hope- 
fully, we can facilitate such a change when we modify the 
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code to allow for the viscosity to depend on the local disk 
parameters as stated above. 

At any heliocentric distance, the timescale for the 
growth of dust grains into planetesimal s is a few thou- 
sand times the local orbital period (jWeidenschillingl 
2000). Therefore, planetesimals formed in inner regions 
of a disk form faster than planetesimals in outer regions. 
Coincidentally, in the inner regions of the planet-forming 
zone our disk in the reference m odel matches with the 
surface density constraints from iDeschl (|2007t ) at early 
times and in the outer planet-forming regions at late 
times. Any planetesimals formed early on in the inner 
region, large enough to decouple from the disk yet not 
so large that they undergo significant migration, will be 
present and available for planet formation at later times. 
These planetesimals would effectively maintain the sur- 
face density of solids at the high levels needed to match 
Desch's (2007) constraints after the gas has been trans- 
ported elsewhere. Farther out in the disk, where plan- 
etesimals form more slowly and at later times, our surface 
density in this region is consistent with surface density 
constraints also at later times. Although we have begun 
the growth of all of the giant planets in our simulations 
at the same time in each model, there is no reason for 
this to be the case. It could very well be that planet 
formation was delayed in one region or another. 

In the simulations LV and LT, we are able to success- 
fully grow the embryos of the outer most giant plan- 
ets within the given ^5—10 Myr time constraint. It 
is important to again stress the uncertainties involved 
in planetary growth from the accumulation of planetes- 
imals. At some point, a transition to oligarchic growth 
and a clearing of the embryos' feeding zones would slow 
embryo growth. On the other hand, other processes ex- 
ist which would increase embryo growth rates. We model 
growth without atmospheres, from a single size distribu- 
tion. We also neglect any tidal effects that could dissipate 
orbital energy from the impacting planetesimals. These 
are but a few of the major uncertainties in core accretion 
models. Any of these and ot hers could e asily alter the 
growth rates by a factor of 2 (jDesch 2003). 

We have also made the simplifying assumption that the 
solids-to-gas ratio is constant with radius. The solids- 
to-gas ratio should increase with radius as exotic ices 
condense in the cold outer regions of the solar nebula. 
The increase of solids beyond the snow line would provide 
more material for planet growth and decrease the growth 
timescale of the giant planets, especially of Neptune and 
Uranus. Photoevaporation of hydrogen and helium from 
the disk would also tend to increase the solids-to-gas ratio 
of the outer disk. Despite these uncertainties, we have 
shown that given reasonable model parameters the cores 
of the outermost giant planets can successfully be built 
in a decreting disk with a truncated outer boundary. It 
is the outward flow of mass and removal at some outer 
radius that provides sufficient mass to the outer regions 
of the disk for planetary core growth. 

We placed our growing embryos in the ad hoc compact 
configuration of the Nice model, but one must keep in 
mind that the Nice model begins after gas has dissipated. 
Planet disk interactions would likely have lead to some 
amount of radial migration while gas was still present 
in sufficient quantities. This would imply that the giant 
planets could have begun to grow elsewhere in the disk 



and migrated to a more compact configuration before the 
gas dissipated. Our placement of the giant planets at the 
locations of the Nice model is most likely incorrect, but it 
is a good proxy for testing planet formation in a compact 
configuration. A more realistic treatment would require 
coupling an A-body code to our viscous evolution code 
and evolving it with migration such that the giant planets 
end up in the compact configuration of the Nice model. 
We plan to investigate this avenue of study in the near 
future. 

It was our aim to simply illustrate how the migra- 
tion timescales are affected by the decreasing surface 
density. One must also keep in mind that type I mi- 
gration is a poorly understood phenomenon. Most of 
the studies to date have investigated type I migration in 
isothermal disks. Accurately coupling an A-body code 
to our viscous disk model with migration will be dif- 
ficult enough without the inherent uncertainties in mi- 
gration processes themselves. It has been shown that 
in non-isothermal disks with high opacities the induced 
net torque m ay have opposite sign and act to push plan- 
ets outward (jPaardckoopcr & Mellcma 2006). That is 
to say, it is uncertain in which direction type I migra- 
tion would force a planet. Furthermore, recent simula- 
tions suggest that the inward scattering of planetesimals 
could drive the outw ard migration of growing embryos 
(|Levison et al.|[20lol ). Type II migration is better under- 
stood. It is likely that the largest giant planet Jupiter 
would quickly fall into the regime of type II migration 
and be swept inward while on the inside of the critical 
radius, where the radial gas flow transitions from inward 
to outward. If, however, the planetary embryo does not 
completely open a gap the type II migration rate ca n be 
reduced or even reversed (jCrida fc Morbidellil 12007') . In 
our model, the critical radius is continually moving in- 
ward. In our reference model, it moves from 38 AU to 14 
AU over the 1.3 Myr of disk evolution. At some point, 
the critical radius should overtake Jupiter and reverse 
the course of its migration outward. It may be that the 
decreasing surface density and outward mass transport 
could save Jupiter from being lost into the Sun. 

9. SUMMARY 

Photoevaporation has for some time now been in- 
voked as a mec hanism for the rapid dispersal of proto- 
planetary disks (|Ballv et al.lll998t I Johnstone et al.lll998l : 
iClarke et all 120011: lAdams et al.1 120041 ). It has recently 
been invoked as a possible mechanism for the truncation 
of the solar nebula in such a fashion as to produce the 
steep surface density profile required to produce the giant 
planets in th e compact configuration of the Nice model 
(jDeschl [20071 ) . We have performed a number of simu- 
lations to test the relative importance of external pho- 
toevaporation versus viscous evolution. We also investi- 
gate whether or not photoevaporative mass loss from the 
outer edge of an evolving protoplanetary disk ca n pro- 
duce t he steep surface density profile posited by IDeschl 
(|2007| ). We find, for reasonable disk parameters, that 
the viscous evolution and photoevaporation play equally 
important roles in determining disk evolution timescales 
and morphology. 

We have adapted a method of solving propagating 
phase change problems and developed a one-dimensional 
viscous disk model that self-consistently tracks the loca- 
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tion of the outer boundary under the influence of pho- 
toevaporation from an external source. Our application 
of the formalism of the Stefan problem to astrophysical 
disks is a novel approach. Our model is, as far as we 
know, the first model to track the location of the outer 
boundary in a fully self-consistent manner. We use radial 
surface density outputs from our viscous disk model to 
model the growth of giant planet embryos in the compact 
configuration set forth in the Nice model. Our model is 
a further ex ploration of t he steady state decretion disk 
proposed by |Desch (2003). With the exception of LV, 
all of our surface density profiles are shallower than the 
surface density profile of Eqn. ([5]), with a typical de- 

pendence on radius of S(r) cx r ' -o.aa. Despite the 
mismatch, our profiles seem conducive to planet forma- 
tio n and do mat ch the surface density constraints derived 
by iDeschl (|2007l ) at certain radii at certain times. 

We present five simulations designed to test the effects 
that varying the strength of the FUV flux and altering 
the strength of the viscosity have on the temporal evolu- 
tion of the disk. It took the disk in LT 3.6 times as long 
to evolve from a mass of 0.07 M Q to a mass of 0.0175 M Q 
as the disk in our reference model, whereas the disk evo- 
lution in HT was faster than of our reference model by 
a factor of 0.45. The evolutionary timescale of LV was 
3.6 times longer than the timescale of our reference model 
and the timescale of HV was 0.81 times shorter than that 
of our reference model. The embryos in LT and LV grow 
faster than the embryos in our reference model, and in 
LV the cores of the three innermost giant planets were 
able to grow to > 10 M ffi within the allocated ^5 — 10 
Myr. 



The strength of the viscosity and the amount of FUV 
flux (envelope temperature) were both able to affect the 
evolutionary timescales of the disks produced in various 
simulations. A small FUV flux or low viscosity were both 
found to produce longer evolutionary timescales than our 
reference model. Both models were more successful than 
our reference model in growing the giant planet cores, 
but they produced very different radial surface density 
profiles. The low FUV flux model produced a radially 
extended disk whereas the low viscosity model produced 
a radially contracted disk with a higher surface density 
in the giant planet forming region. The differences in 
these two models may provide a natural explanation for 
the location of the outer edge of the solar system. 

The disks produced with our numerical model are all 
characterized by outward mass transport, mass loss at 
the outer edge, and a truncated outer boundary. The 
outer boundary is characterized by substantial mass loss 
due to photoevaporative heating. This mass loss drives 
outward mass flow from the critical radius to the outer 
edge of the disk. The transport of mass from small radii 
to large can potentially prevent the rapid inward migra- 
tion of Jupiter and Saturn, while at the same time supply 
enough mass to the outer regions of the disk for the for- 
mation of Uranus and Neptune. 

The authors thank John Bally for his illuminating dis- 
cussions and insightful comments. We thank the anony- 
mous referee for helpful comments that improved the 
quality of the paper. This work was supported by the 
Cassini Project. 
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VARIABLE SPACE GRID METHOD 

A well-posed problem in the material sciences, called the Stefan Problem, deals with propagating phase changes 
typically considered in the context of melting/freezing and heat ablation (Ozisik .1980). We have developed a 1-D 
model of a protoplanetary disk that includes viscous diffusion and photoevaporation at the outer boundary using 
numerical techniques developed to solve the Stefan Problem. By adapting the Stefan problem to astrophysical disks, 
our 1-D numerical model self-consistcntly tracks the location of the outer boundary. This is a novel approach to 
modeling the solar nebula. 

To solve t he Stefan problem, iKutluav et all (|1997l ) adopt a numerical method with a variable space grid (VSG) first 
proposed bv lMurray &: Landisl (|195^ ). The VSG method employs a fixed number of grid points with a variable grid 
size at each time step. This method involves solving two coupled differential equations at each time step, one for the 
location of the outer boundary and one for the diffusive evolution of the disk. Once the location of the outer boundary 
is found the abscissa is rescaled and the diffusive evolution calculated. 

We will now briefly describe the VSG method. For a non-dimensional diffusion equation, 



with the boundary condition 
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where s(t) is the location of the boundary at time t. Eqn. (| Al|) can be differentiated with respect to time and, for the 
zth grid point, 
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By substituting Eqn. (|A4I) into Eqn. (|A3[) . the heat equation (Eqn. (|A1[1 ) can be reformulated as 

m_j^ds(f)_dU d 2 U 
dt s(t) dt dx dx 2 

Various methods, such as similarity solutions, have been used to solve Stefan problems analytically. We used these 
solutions to test our numerical cod e. One such test was done on the following system, a Stefan problem of transient 
heat conduction in a melting slab (jBluman fc Coldll974[ ). The one-dimensional, finite slab is insulated at x = and 
has a propagating phase change at x = X(t), where heat flows into the melting face at a rate H(t). 

du d 2 u , , , . . 

Tt=^ Q<x<x V ( ) 

with the boundary condition 

Given the following boundary and initial conditions 

u(0,t) = 0, x = (A8) 



u(X(t),t) = 0, x = X(t) (A9) 

and 

u(x,0) = g(x) (A10) 
and assuming X(t) = 1 — t, the exact solution is of a self-similar form 

u(x, t) = exp(7r ) /2 cxpl — J (All) 

where £ = x/X(t). 

The exact analytic, self-similar solution was used to check the numerical VSG method. We tested our code for 
convergence against the exact analytical solutions by increasing the spatial grid resolution. These tests were all done 
with the same sized time steps. The results of these tests have been tabulated in Table [2] The numerical results are 
in good agreement with the exact solution and exhibits the expected convergence as the number of grid point, N, 
increases. 



TABLE 2 

Analytic test for the VSG method. Results of the tests for convergence against exact, analytic self-similar solution 

AT A FINAL TIME OF t = 0.25 ARE SHOWN. 



x/X(t) 


u (actual) 


u(N= 50) 


u (N = 100) 


u(N = 200) 


0.0 


0.00000000 


0.00000000 


0.00000000 


0.00000000 


0.2 


0.02547849 


0.02540785 


0.02546055 


0.02547405 


0.4 


0.04216313 


0.04204546 


0.04213325 


0.04215575 


0.6 


0.04377427 


0.04365066 


0.04374288 


0.04376651 


0.8 


0.02851226 


0.02843037 


0.02849147 


0.02850713 


1.0 


0.00000000 


0.00000000 


0.00000000 


0.00000000 



As a further test of convergence on the problem at hand we have completed a number of simulations with a variety 
of grid sizes. We have checked for both spatial and temporal convergence. Figure [TT] shows the spatial convergence as 
the number of grid spaces increases. The convergence has been calculated by differencing the surface density of each 
of the lower resolution simulations from the surface density of the highest resolution simulation (N — 800) and then 
normalizing by the innermost available grid space. The computed convergence is shown with plus symbols connected 
by solid lines. For comparison, the expected 1/N 2 convergence is shown over-plotted with x's connected by dotted 
lines. It can clearly be seen that the actual convergence is very close to the expected convergence. We have chosen 
to use N — 200 for our number of grid spaces. This allows for simulations that complete in a reasonable amount 
of time and is acceptably accurate, to within less than 0.5% of the highest resolution simulation. Due to the large 
uncertainties in many model parameters, we feel this level of convergence is acceptable. 

We have also investigated the temporal convergence of our model. For this test we have divided the end times of 
each simulation by the end time of the highest resolution simulation. Again, the highest resolution simulation has 
N = 800. The temporal convergence was calculated by dividing the final time of each lower resolution simulation by 
the final time of the highest resolution simulation. It can be seen in Figure [12] that the final time of the simulation 
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Fig. 11. — The spatial convergence of a number of simulations with various numbers of grid spacings. The calculated convergence is 
shown with plus symbols connected by solid lines. The expected 1/N 2 convergence is shown with x's connected by dotted lines. In all of 
our simulations we use N = 200 grid spaces. 




200 400 600 800 

grid resolution [N] 

Fig. 12. — The temporal convergence of a number of simulations with various numbers of grid spacings. In all of our simulations we use 
N = 200 grid spaces. 



with N = 200 is within 2% of the final time of the highest resolution simulation. As with the spatial convergence, we 
feel that this is sufficient considering the large uncertainties in many of our model parameters. 
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